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In this letter, a generalization of pairwise models to non-Markovian epidemics on networks is 
presented. For the case of infectious periods of hxed length, the resulting pairwise model is a 
system of delay differential equations (DDE), which shows excellent agreement with results based 
on explicit stochastic simulations of non-Markovian epidemics on networks. Furthermore, we 
analytically compute a new 7?.o-like threshold quantity and an implicit analytical relation between 
this and the final epidemic size. In addition we show that the pairwise model and the analytic 
calculations can be generalized in terms of integro-differential equations to any distribution of the 
infectious period, and we illustrate this by presenting a closed form expression for the final epidemic 
size. By showing the rigorous mathematical link between non-Markovian network epidemics and 
pairwise DDEs, we provide the framework for a deeper and more rigorous understanding of the im¬ 
pact of non-Markovian dynamics with explicit results for final epidemic size and threshold quantities. 


Networks have provided a step change in modelling 
complex systems ranging from man-made to natural ones 
[1-3] . The study of disease transmission on networks has 
particularly benefitted from this modelling paradigm by 
uncovering the role and impact of contact heterogeneity 
and clustering [3], to name just a few. While networks 
provide a clear departure from classic compartmental 
models, the role of mean-field models remains crucial. 
These offer a reliable way to obtaining analytical res¬ 
ults and thus uncovering the interplay between network 
properties and the dynamic processes on networks. For 
example, the epidemic threshold [4, 5] and final epidemic 
size [6] can be given in terms of explicit or implicit math¬ 
ematical expressions which clearly illustrate how network 
and disease parameters combine. 

Probably the most widely spread and well-known 
mean-field model for network epidemics is the degree- 
based mean-field (DBMF) model, also known as hetero¬ 
geneous mean-field [3, 4]. Similarly, pairwise models [6- 
9] continue to provide a fruitful framework for model¬ 
ling dynamic or adaptive networks involving epidemics 
[8, 10], social interactions [11] and ecological systems [9]. 
Such models come with the added benefit of some degree 
of analytical tractability and the means toward expli¬ 
cit analytical quantities such as the basic reproduction 
number and final epidemic size [6]. Recently, however 
there is renewed interest in modelling non-Markovian 
processes, such as epidemics on networks [12-15], ran¬ 
dom walks [16] and temporal networks [17]. This recent 
burst of research focusing on non-Markovian dynamics 
is strongly motivated by empirical observations. These 
show that for many real world settings, the Markovian 
framework is not satisfactory in describing temporal stat¬ 


istics, such as time intervals between discrete, consecut¬ 
ive events. Examples include inter-order and inter-trade 
durations in financial markets [18], socio-networks [19], 
or individual-to-individual contacts being dynamic [17]. 
In the context of epidemiology, the period of infectious¬ 
ness has paramount importance [20, 21], and researchers 
departed from the simplifying assumption of exponential 
distributions by approximating the empirical distribu¬ 
tion of observed infectious periods of various diseases by 
log-normal and gamma (smallpox [22, 23]), fixed-length 
(measles [24]) or Weibull distributions (ebola [25]). The 
reliable tools and mathematical machinery of Markovian 
theory do not translate directly to modelling and analysis 
of non-Markovian systems, and this is the main source of 
many challenges. 

In this letter, we present the first analog of pairwise 
models for non-Markovian epidemics, and show that this 
is equivalent to a set of delay differential equations which 
(a) show excellent agreement with simulation and (b) 
allows us to give an implicit analytic expression for the 
final epidemic size, as well as to define a new TT-g-like 
quantity which emerges naturally from our calculations. 

We consider an undirected and unweighted network 
with N nodes and an average degree n. Each node can 
be susceptible (S), infected (J) and recovered (R). For 
Markovian epidemics, with transmission rate r and re¬ 
covery rate 7, the epidemic is well approximated by the 
pairwise model [6] given below 

[^] = -t[^/], [/] = r[5/] - 7[/], [SS] = -2t[SSI], 

[5/] = t[SSI] - t[ISI] - t[5'/] - 7 [S'/], 

where [X], [XY] and [XYZ] are the expected num¬ 
ber of nodes in state X, links in state X — Y and 
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triples in state X — Y — Z, respectively. Considering 
the network at a given time, then counting amounts to 
W = EhXi, [xy] = and [XYZ] = 

k=i^^YjZkg^jgjk, where X,Y,Z G {S,I,R}, and 
G = is the adjacency matrix of the net¬ 

work such that gu = 0, gij = gji and gij = gji = 
1 if nodes i and j are connected and zero otherwise. 
Moreover, Xi returns one if node i is in state X and zero 
otherwise. The dependence on higher order moments can 
be broken by using that [XSY] = [6]- Ap¬ 

plying this leads to the following self-consistent system 


[5] =[/] = r[,5/] - 7[/], 

n- 1 [SS][SI] 


[SS”] = -2r- 


[^] 


^ - 1 [SI][SI] 


n [S'] 
-r[S/]- 7 [S/]. 


[* 5 ] 


( 1 ) 


By applying the closure at the level of pairs, [XY] = 
n[X]^, system (1) reduces to the classic compartmental 
SIR model. 


• T) • n 

S = -r-S/,/ = r-S/-7/. 


( 2 ) 


We wish to apply the previous approach to the case 
when the recovery time is not exponentially distributed. 
First, a fixed infectious period, denoted by cr, is con¬ 
sidered, and the derivation of the pairwise model from 
first principles is illustrated. We show that the non- 
Markovian dynamics can be described by a delay differ¬ 
ential equation with constant delay. The infection pro¬ 
cess is assumed to be Markovian, thus the equation for 
[S'] is the same as before, namely [S]{t) = —r[S/](t). 
The number of infected nodes at time t is replenished by 
T[S/](t) and is depleted by r[S/](t — cr), and this yields 
[I]{t) = T[SI](t) — T[SI]{t — a). The equation for the 
number of S — S links is the same because the infection 
process is Markovian, see (1). In a similar manner, the 
number of S — / links is replenished by > 

which is the rate of depletion of S— S links. Furthermore, 
depletion occurs due to the infection within S — I pairs, 
T[S/](t), and due to the infection of a central S node in 
an / — S — / triple, [SI](t) [SJ](t). On the other 

hand, there are S — I links, which survive the time in¬ 
terval cr, that will be removed due to the recovery of the 
I node. However, one needs to account for the removal 
of S' — / links which were created precisely a times ago. 
Naively, one would believe that this term is simply pro¬ 
portional to ■ However, one must 

account for the fact that in the time interval (t — a, t) an 
S — I link could have been destroyed either due to within 
pair infection or by infection of the S node from outside. 
Hence, it is obvious that a discount factor needs to be 
determined to account for this effect. To calculate this 


factor, S — I links, that are created at the same time, 
are considered as a cohort denoted by x, and we model 
infection within and from outside by writing down the 
following evolution equation, 

x'{t) = -^^^^j^[SI]{t)x{t)-Tx{t), (3) 

where, the first term denotes the ‘outer’ infection of the S 
node, while the second term stands for ‘inner’ infection of 
the S node. We note that the outside infection is simply 
proportional to the probability that an S node with an 
already engaged link has a further susceptible neighbour, 
solution of equation (3) in time interval 

[t — a, t] is 

x{t) = x{t - a)e- fL.{^^lSi]iu)+T)du^ 

and this provides the depletion or discount rate of S' — / 
links. In this case, x{t — a) = ^ 

which is the replenishment of S — / links. Therefore, 
summarising all the above, the pairwise DDE for the non- 
Markovian case is 


[S]it) = -T[SI]it), m) = r[S/](t) - r[SI]{t - a) 
[SS](t) = [Sl]{t) = -T[S/](t) 


[sm 


T{n — I) 


[simisim 


n-l [SS]{t)[SI]{t) 


(4) 


n[S](t) n [S](t) 

_ _ n- 1 [SS](t- g)[S/](f- cr) _ ft 
n [S](t —cr) 


This system is now the main subject of our investigation 
from analytical and numerical point of view. Similarly 
to Markovian case, the non-Markovian mean-field model 
for fixed infectious period is 

• T) 

S{t) = -r-S{t)I{t), 

Ht) = T-^S{t)Iit) - r^S{t - a)I{t - cr). (5) 

The most important qualitative results for SIR mod¬ 
els are the explicit formula of basic reproduction number 
and an implicit equation for the final epidemic size. In 
what follows, we introduce a general concept for the re¬ 
production number associated to pairwise models, and we 
refer to this as the pairwise reproduction number. Using 
this concept, the final size relations for the above mean- 
field, classic pairwise and DDE-based pairwise models 
are derived. Reproduction numbers play a crucial role 
in mathematical epidemiology, so we begin by investig¬ 
ating these. The basic reproduction number Rq denotes 
the expected number of secondary infections caused by a 
‘typical’ infected individual during its infectious period 
when placed in a fully susceptible population, which is a 
definition understood at the level of nodes (individuals). 
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On the other hand, the pairwise model is written at the 
level of links and describes the dynamics of susceptible 
{S — S) and infected {S — I) links. This fact gives us an 
opportunity to define a new type of reproduction num¬ 
bers, which we call pairwise reproduction number and 
denote it by TZq. More precisely, we distinguish the fol¬ 
lowing two useful quantities: (a) the basic reproduction 
number is the expected lifetime of an / node multiplied 
by the number of newly infected nodes per unit time, 
and (b) the pairwise reproduction number is the expec¬ 
ted lifetime of an S' — / link multiplied by the number of 
newly generated S — I links per unit time. 

An infected node is removed due to its recovery, thus 
in general the expected lifetime is the expected value of 
a random variable X corresponding to the distribution 
of the length of infectious periods. In contrast, an S — / 
link can be removed due to the recovery of the / node but 
also due to the infection of the S node. Therefore, the 
expected lifetime of the S — I link is the expected value 
of the minimum of two random variables. If we assume 
that the process of infection along such a link has density 
function fi with survival function and the process of 
recovery has density function fr with survival function 
then, denoting by Z the random variable dehned by 
the lifetime of an S' — / link, we have 

/ CO 

X + fr{x)^i{x)) dx. (6) 

-OO 

From the assumption that the infection time along S — 
I links is exponentially distributed, the number of newly 
infected nodes per unit time is ;^t[S]o in the mean-field 
model, and the expected number of newly infected links 
is = T^^[S]o in the pairwise model, where 

we used the approximation [SS]o = ;^[S]g. 

We illustrate how to use the formula (6) in the case 
of fixed length infectious period (cr). In this case, the 
survival function is 


Ut) 


1 if 0 < t < cr, 
0 if t > cr. 


and the density function fr{t) is the Dirac-delta 6{t — a). 
Using fundamental properties of the delta function, we 
have 


/ OO pCO 

xfi{x)^r{x)dx + / xfr{x)^i{x)da 

-OO J —OO 


1 — e 


= —ere 


-|- ere 



TZo 

K 

Markovian 

n T Q 
iV 7*^0 

n — 1 T rci 

JV T+7l'^lo 

Fixed 


^(l-e-"'^)[5]o 

General 

frE(X)5o 

V(i-^[A](^))[5]o 


Table I. Basic and pairwise reproduction numbers for different 
recovery distributions. C[fr\{r) denotes the Laplace trans¬ 
form of fr, the density of the recovery process, at r. 


This provides a very general result, which in many cases 
leads to an explicit analytical result for TZq, see Table I. 

For the standard Markovian mean-held model, the pro¬ 
cess of calculating the hnal epidemic size is well-known. 
From Eq. (2), we evaluate dl/dS and integrate it to 
obtain 



Using that TZq = ^j^Sq, we have 



The hnal epidemic size (i.e. the total number of infec¬ 
tions) can be easily computed by using Roo = N—Soo- In 
the non-Markovian case, the calculations (which are in¬ 
cluded in the supplemental material) are rather different 
and the resulting hnal size relation is 


In 



n „ 



( 8 ) 


As in this case TZq = r^aSo, the hnal size relation (8) 
shows the ‘standard’ form of (7). The dynamical sys¬ 
tems (1) and (4) can be manipulated conveniently to de¬ 
rive an analytic relation between the hnal epidemic size 
and the basic reproduction number. This is known for 
the Markovian case but it is a new result for the non- 
Markovian one. While the full derivation for the non- 
Markovian case is given in the supplemental material, the 
main steps of the calculations are: (a) hnd an invariant 
to reduce the dimensionality of the system, (b) integrate 
the equation for [5'/](t), (c) integrate the equation for 
[S]{t) on [0,oo) and (d) employ algebraic manipulations 
to obtain the hnal size relation. Following this procedure 
yields 

^^ = ^(l-e-"")[5]o(s^-l), (9) 

n—1 


and multiplying this result by the number of newly gen¬ 
erated S — I links, the formula in Table I for TZq fol¬ 
lows. More importantly, it is noteworthy to highlight 
the general result that IE(2'), upon using that the in¬ 
fection process is Markovian, reduces to evaluating the 
Laplace transform of the density of the recovery time. 


where Soo = and the attack rate is simply 1 — Soo- 
Using the same technique for the Markovian case leads 
to 

^^[s|„(,jr-i). (10) 


1 

n — 1 
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Upon inspecting the two relations above, the following 
important observations can be made. First, the implicit 
relation between final size and TZ^ is conserved between 
the Markovian and non-Markovian model. Moreover, 
upon using the values of TZq as given in Table I, equations 
(9) and (10) can be cast in the following general form 

(11) 

n—1 

In fact we conjecture that this relation will hold true for 
pairwise models with different infectivity period profiles. 
The second observation is that taking the limit of n —>■ oo 
in (11) gives rise to 

In(s^) = 7^g(so, - 1), (12) 


which is equivalent to the ‘standard’ form of (7) 

To test the validity of our model we implemented an 
algorithm to simulate the non-Markovian SIR process 
with arbitrary recovery times, and considered random 
networks with N = 1000 nodes. In Fig. l(a,b) homogen¬ 
ous and Erdds-Renyi random networks are considered, 
respectively. Here, the mean of 100 simulations is com¬ 
pared to the solution of system (4). The agreement is 
excellent for homogenous networks even for low degrees. 
Despite the pairwise model not accounting explicitly for 
the network’s degree distribution, the agreement is sur¬ 
prisingly good for relatively dense Erdds-Renyi networks. 

In Fig. 1(c) we compare and contrast the differences 
between simulations, mean-field and pairwise models for 
the non-Markovian case. For denser networks {{k) = 
15), both models perform well with the pairwise yield¬ 
ing a better agreement with output from simulation. 
However, the difference is striking for sparser networks 
((fc) = 5), where the mean-field approximation performs 
very poorly, while the pairwise DDE model leads to good 
agreement even in this case. 

In Fig. 1 (d), analytic final size relations are tested 
against simulation results for a range of different infec¬ 
tious period distributions, all sharing the same mean. 
The horizontal axis corresponds to the transmission rate 
T, and the plot highlights the threshold dynamics, as well 
as the necessity to correctly model the recovery time dis¬ 
tribution in order to avoid under or over estimation of 
the final epidemic size. Based on Table I, the analytical 
expressions for TZq can be computed for the Markovian 
(or exponential), fixed and gamma-distributed recovery 
times. These values are 


'^o,r(U2) ^ 

"^o,r(2,i) = c ^1 - 
where c = 


1 


VTT^ 

4 




0,Exp(l) 


= C 


T + l 


^0,Fixed(l) “ ^ ® ) ’ 


(2 + r)2 

, and satisfy the following inequality 


nl^,. (13) 


o.r(U2) — o,Exp(i) — '''o,r(2,i) 


We note that (a) all recovery time distributions have the 
same mean 1 and (b) the variances satisfy the converse 
inequality, with higher variance in recovery time (i.e. 2, 
1, 1/2 and 0) giving a smaller value, despite t be¬ 
ing fixed. The overall agreement between the analytic 
results of the pairwise model and the stochastic simula¬ 
tions is excellent and confirms the validity of our final 
size relations. The inset in Fig. 1 (d) illustrates how the 
final epidemic size depends on the pairwise reproduction 
number, and shows that the same value of TZ^ produces 
the same attack rate, regardless of the distribution from 
where it is originated from, in accordance with our for¬ 
mula (11). 

We have introduced a generalization of pairwise models 
to non-Markovian epidemics with fixed infectious period. 
The resulting model is a system of delay differential equa¬ 
tions with constant delay and we have provided as a 
full as possible analytical and numerical analysis of this 
model and benchmarked its performance against expli¬ 
cit stochastic network simulations. We have presented 
a new concept of reproduction numbers introducing the 
pairwise reproduction number TZq and have derived the 
final epidemic size relation for non-Markovian mean-field 
and pairwise DDE models. 

The numerical solution of the non-Markovian pairwise 
DDE shows excellent agreement with results based on 
explicit stochastic network simulations and sheds some 
light on the impact of non-Markovianity. More import¬ 
antly, via the analytic results we can gain insight how 
and where non-Markovianity enters and impacts upon 
important epidemic descriptors. 

The model and results in this paper should provide 
a framework for deeper and more comprehensive ana¬ 
lysis of non-Markovian processes on networks and these 
should not be restricted to epidemics with fixed delays. 
Preliminary investigations indicate that our model can 
be extended to consider arbitrary recovery time distribu¬ 
tions. In this case, the resulting model is a more complex 
integro-differential equation requiring a more challenging 
and elaborate analysis. Nevertheless, it turns out that 
the hnal epidemic size relation, upon assuming a general 
probability distribution with density function (fr), yields 

"-^='^il-C[fr]ir))[S]o{sT-l), (14) 

n—1 

which agrees with the general equation suggested in (11). 
For recovery of fixed length, relation (14) reduces to 
(9). The validity of our general final size relation (14) 
was tested also for different gamma distributions, see 
Fig. 1(d), showing a strong predictive power for general 
non-Markovian epidemics on networks. The difficulty of 
modelling non-Markovian processes is well known, but 
our current framework can pave the way for identifying 
fruitful links between different areas of delay differential 
equations, stochastic processes, dynamical systems and 
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Figure 1. Simulations of non-Markovian epidemics on networks with N — 1000 nodes: (a) solid lines show the numerical 
solution of (4) and the circles/squares/diamonds correspond to simulations for homogeneous networks with (k) = 5/10/15, 
respectively; (b) the same as before but for Erdos-Renyi random networks with (k) = 5/10/15; (c) the solid and dashed lines 
show the numerical solution of pairwise (4) and mean-field (5) models, respectively and, for homogeneous networks with (k) = 5 
and (k) = 15. For (a), (b) and (c) the transmission rate is r = 0.55 and the infectious period is fixed, ct = 1. Finally, (d) 
the diamonds/circles/squares correspond to numerical simulations using homogeneous network with (k) — 15 and using fixed 
and two different but gamma distributed infectious periods (o - shape a = 2, scale /?=!,□- shape a = |, scale /3 = 2), 
respectively. The solid lines correspond to the analytical final epidemic size for fixed (9) and general (14) infectious periods. 
The inset shows the analytical and the simulated final epidemic sizes plotted against the pairwise reproduction number. 


epidemiological models. 
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